************************************************************
****************** REPLICATION CODE ************************
************************************************************
* For: Rigterink, The Wane of Command, APSR				   *
* Simulation of reporting bias							   *
************************************************************
* Content of file:										   *
/*														   *
1. Preliminaries and setting directory					   *
2. Finding distribution of terrorist attacks after hit	   *
3. Main result											   *
4. Setting parameters for each terrorist group			   *
5. Simulation of reporting bias							   *
6. Outputting results table								  */
************************************************************

/*----------------------------------------------------------
---------- Preliminaries and setting directory -------------
----------------------------------------------------------*/

clear all
version 14
set more off

*******************************
*** Set your directory here ***
*******************************
cd "[directory]"

cap log close
log using log\Rigterink_drones_replication_appendix_reportingbias.log, replace

use "dta/Rigterink_drones_replication_maindataset.dta", clear


/*----------------------------------------------------------
-- Finding distribution of terrorist attacks after a hit ---
----------------------------------------------------------*/

* Indicators for periods six months after a hit / miss
gen posthit = 1 if bigfishdied==1 | L1.bigfishdied==1 | L2.bigfishdied==1 | L3.bigfishdied==1 | L4.bigfishdied==1 | L5.bigfishdied==1 | L6.bigfishdied==1
gen postmiss = 1 if (bigfishtarget==1 & bigfishdied==0) | (L1.bigfishtarget==1 & L1.bigfishdied==0) | (L2.bigfishtarget==1 & L2.bigfishdied==0) | (L3.bigfishtarget==1 & L3.bigfishdied==0) | (L4.bigfishtarget==1 & L4.bigfishdied==0) | (L5.bigfishtarget==1 & L5.bigfishdied==0) | (L6.bigfishtarget==1 & L6.bigfishdied==0)
replace postmiss = . if posthit==1

* For each group, test poisson, nbreg and zinb for best fit
forvalues i = 1(1)13 {

	qui sum bigfishdied if groupid==`i'
	scalar anyhits = r(max)
	
	if anyhits==1 {
	
		cap poisson terratt if groupid==`i' & posthit==1, iterate(250)
		scalar ll_pois`i' = e(ll)

		cap nbreg terratt if groupid==`i' & posthit==1, iterate(250)
		scalar ll_nb`i' = e(ll)
		
		cap zinb terratt if groupid==`i' & posthit==1, inflate(_cons) , iterate(250)
		scalar ll_zinb`i' = e(ll)
		
		di "Group number `i'"
		
		di "NB vs ZINB"
		di "LL NB `i': " ll_nb`i'
		di "LL ZINB `i': " ll_zinb`i'
		di "Prob > chi2 = " chi2tail(2, 2*(abs(ll_zinb`i'-ll_nb`i')))
	
		di "NB vs POISSON"
		di "LL NB `i': " ll_nb`i'
		di "LL POISSON `i': " ll_pois`i'
		di "Prob > chi2 = " chi2tail(2, 2*(abs(ll_pois`i'-ll_nb`i')))
		
		
	}
}

* For three out of 5 groups with a hit, evidence that NB outperforms Poisson
* Evidence for no groups that ZINB outperforms NB


/*----------------------------------------------------------
---------------------- Main result ------------------------
----------------------------------------------------------*/

* Producing main results in main article
reghdfe terratt L(0/6).bigfishdied L(0/6).bigfishtarget L(0/6).dronestrike F(1/6).bigfishtarget F(1/6).bigfishdied F(1/6).dronestrike, absorb(i.groupid i.periodid) vce(, bw(12))
matrix actual = r(table)

* Determining if any of the coefficients are significant at the 5% level
forvalues i = 2(1)7 {
	local coef`i' = actual[1,`i']
	local pval`i' = actual[4,`i']
}

 if (`coef2'>0 & `pval2'<0.05) | /// 
(`coef3'>0 & `pval3'<0.05) | ///
(`coef4'>0 & `pval4'<0.05) | ///
(`coef5'>0 & `pval5'<0.05) | ///
(`coef6'>0 & `pval6'<0.05) | ///
(`coef7'>0 & `pval7'<0.05) {
	local significant = 1
}


/*----------------------------------------------------------
------ Setting parameters for each terrorist group ---------
----------------------------------------------------------*/

* Matrix with parameters
matrix para = J(13, 5, .)

* For each group, determine parameters distribution attacks using nbreg
forvalues i=1(1)13 {

	qui sum bigfishdied if groupid==`i'
	scalar anyhits = r(max)
	
	if anyhits==1 {
		qui nbreg terratt_unlog if groupid==`i' & posthit==1, iterate(100)
	}
	
	if anyhits==0 {
		qui nbreg terratt_unlog if groupid==`i', iterate(100)
	}
	
	matrix para[`i', 1]= e(b)
	matrix para[`i', 3]= e(alpha)
	matrix para[`i', 4]= 1/para[`i', 3]
	matrix para[`i', 5]= exp(para[`i', 1])
	
}


/*----------------------------------------------------------
-------------- Simulation of reporting bias ----------------
----------------------------------------------------------*/

set more off
set matsize 11000

forvalues R=1(1)100 {
	
	local count = 1
	
	di "Round: " `R'
	timer on 1
	
	* Setting seed
	local seed = 19860223+`R'
	set seed `seed'
	
	* Generating simulated dependent variable
	qui gen simdep=.
	
	forvalues i = 1(1)13 {
		qui gen xg = rgamma(para[`i', 4], para[`i', 3]) if groupid==`i'
		qui gen xbg = xg*para[`i', 5] if groupid==`i'
		qui replace simdep = rpoisson(xbg) if groupid==`i'
	
		qui drop xg xbg
	}
	
	qui replace simdep=0 if simdep==.
	
	*  Repetitions for each simulated dataset
	forvalues N=1(1)10 {
	
		local seed = `seed'+`N'
		set seed `seed'
	
		di "Sub-round: " `N'
		
		* Probability of noticing terrorist attack after miss
		forvalues Pmiss = 0.1(0.05)1 {
			
			* Probability of noticing terrorist attack in absence of leader-involved drone strike
			forvalues Pnone = 0.1(0.1)1 {
				
				* Assuming attacks are more noticeable after a drone strike then after no drone strike
				if `Pmiss'>`Pnone' {
				
					* Simulating whether attack gets noticed
					qui gen noticed = simdep if posthit==1
					qui replace noticed = rbinomial(simdep, `Pmiss') if postmiss==1
					qui replace noticed = rbinomial(simdep, `Pnone') if posthit==. & postmiss==.
					qui replace noticed = 0 if simdep==0
					qui replace noticed = ln(noticed+1)
			
					qui reghdfe noticed L(0/6).bigfishdied L(0/6).bigfishtarget L(0/6).dronestrike F(1/6).bigfishtarget F(1/6).bigfishdied F(1/6).dronestrike, absorb(i.groupid i.periodid) vce(, bw(12))
					matrix simulated = r(table)
					qui testparm L1.bigfishdied L2.bigfishdied L3.bigfishdied L4.bigfishdied L5.bigfishdied L6.bigfishdied
					local simpval = r(p)
					
					* Determining if any of the coefficients are significant at the 5% level
					forvalues j = 2(1)7 {
						local coef`j' = simulated[1,`j']
						local pval`j' = simulated[4,`j']
					}

					if (`coef2'>0 & `pval2'<0.05) | /// 
					(`coef3'>0 & `pval3'<0.05) | ///
					(`coef4'>0 & `pval4'<0.05) | ///
					(`coef5'>0 & `pval5'<0.05) | ///
					(`coef6'>0 & `pval6'<0.05) | ///
					(`coef7'>0 & `pval7'<0.05) {
					local significant = 1
					}
					else {
					local significant = 0
					}
					
					* Storing result in matrix
					matrix result = (`R', `N', `Pmiss', `Pnone', `significant', `simpval', `coef2', `coef3', `coef4', `coef5', `coef6', `coef7', `pval2', `pval3', `pval4', `pval5', `pval6', `pval7')
					
					if `count'==1 {
						matrix allresults`R' = result
					}
					else {
						matrix allresults`R' = (allresults`R' \ result)
					}
				
					drop noticed
					local count = `count'+1 
				}
			}
		}
	}
			
	drop simdep
	timer off 1
	di "Time elapsed for Round:"
	timer list 1
}
	

* Storing all results in matrix and converting to dataset

set more off

forvalues R=1(1)100 {
	clear
	svmat allresults`R'
	rename (allresults`R'1 allresults`R'2 allresults`R'3 allresults`R'4 allresults`R'5 allresults`R'6 allresults`R'7 allresults`R'8 allresults`R'9 allresults`R'10 allresults`R'11 allresults`R'12 allresults`R'13 allresults`R'14 allresults`R'15 allresults`R'16 allresults`R'17 allresults`R'18) (round subround p_miss p_none significant Fstat coef1 coef2 coef3 coef4 coef5 coef6 pval1 pval2 pval3 pval4 pval5 pval6)

	if `R'!=1 {
		append using dta/Rigterink_drones_replication_reportingbias.dta
	}

	save dta/Rigterink_drones_replication_reportingbias.dta, replace
	
}

save dta/Rigterink_drones_replication_reportingbias.dta, replace


/*----------------------------------------------------------
---------------- Outputting results table ------------------
----------------------------------------------------------*/
	
use dta/Rigterink_drones_replication_reportingbias.dta, clear

* Collapsing by 5th and 95th percentile
for var coef?: gen X_p5 = X \ rename X X_p95
collapse (mean) significant (p5) coef?_p5 (p95) coef?_p95, by(p_miss p_none)
drop *_p*

* Rounding
replace p_none = round(p_none*10, 1)

* Reshaping
reshape wide significant, i(p_miss) j(p_none)

* Extra observation for heading
set obs `=_N+1'
replace p_miss = 100 if p_miss==.

* Putting observations in order
gsort - p_miss
order p_miss significant9 significant8 significant7 significant6 significant5 significant4 significant4 significant3 significant2 significant1

* Rounding
replace p_miss = round(p_miss*100, 1)
tostring p_miss, replace

* Making header
replace p_miss="P(report | miss)" if p_miss=="10000"

* Making left column
forvalues i=1(1)9 {
	label var significant`i' "0.`i'"
}

label var p_miss "P(report | none)"

gen comma = "0."
egen p_miss_temp = concat(comma p_miss) if _n!=1
replace p_miss = p_miss_temp if _n!=1
drop comma p_miss_temp

* Outputting table
texsave using "tables/Rigterink_ATabB1.tex", replace nofix hlines(1) varlab frag title("Simulation of reporting bias") location(h) marker("Atab:reportbias") footnote("This table displays the share of simulated regressions with at least one coefficient statistically significant at the 5\% level at the given probability that media report an attack by a terrorist group in the six months after a drone miss on its leader, and given no drone attempt on its leader respectively. Probability of media reporting an attack by a terrorist group in the six months after a drone hit on its leader is set to 1.")

*******************
*** END OF FILE ***
*******************
